Python for samfunnsøkonomi
Forelesningsnotater - Espen Sirnes
10 - statsmodels
Statsmodels er den mest brukte pakken for å kjøre statistiske analyser i Python. I denne forelesningen skal vi kun se på den enkleste typen regresjon - Minste Kvadraters Metode (Ordinary Least Squares) - eller OLS som vi vanligvis kaller det.
For å bruke statsmodels må vi ha et datasett. Vi skal bruke et datasett fra gapminder.com. Datasettet viser sammenhengen mellom BNP og forventet levealder på en illustrerende måte.
Vi begynner med å laste inn datasettet fra https://titlon.uit.no/hht/data/gapminder.csv:
import pandas as pd
import requests
# URL of the CSV file
url = "https://titlon.uit.no/hht/data/gapminder.csv"
#sometimes this doesn't work on mac
#g = pd.read_csv(url)#reading data
# Local file path to save the downloaded CSV
local_file_path = "gapminder.csv"
# Download the CSV file using requests
response = requests.get(url)
with open(local_file_path, 'wb') as f:
f.write(response.content)
# Read the locally saved CSV file into a Pandas DataFrame
g = pd.read_csv(local_file_path)
g
Vi plotter så forventet antall leveår mot produksjon (BNP) per innbyger, slik vi har lært i 3 -matplotlib:
from matplotlib import pyplot as plt
fig,ax=plt.subplots()
#adding axis lables:
ax.set_ylabel('Forventet antall leveår')
ax.set_xlabel('BNP per innbygger')
#plotting the function:
ax.scatter(g['gdp_cap'], g['life_exp'], label='Observasjoner')
ax.legend(loc='lower right',frameon=False)
Grafisk ser det ut til at forventet levealder øker mye for lavt BNP per innbygger. Det tyder på at det er den prosentvise økningen i BNP som er avgjørende. Altså at en økning i inntekt fra for eksempel hunder kroner dagen til to hundre kroner, har en mye større effekt enn en økning fra kr 2 000 til kr 2 100. Ved å konvertere bnp til log blir den prosentvise økningen konstant langs x-aksen:
import numpy as np
from matplotlib import pyplot as plt
fig,ax=plt.subplots()
#adding axis lables:
ax.set_ylabel('Forventet antall leveår')
ax.set_xlabel('BNP per innbygger')
#plotting the function:
ax.scatter(np.log(g['gdp_cap']), g['life_exp'], label='Observasjoner')
ax.legend(loc='lower right',frameon=False)
Vi ser at vi nå får en penere spredning, og det ser ut til å være en sammenheng. Vi skal nå forsøke å tegne en "regesjonslinje" blant prikkene som er slik at avstanden til prikkene er så liten som mulig, i gjennomsnitt. Det er denne metoden som kalles "Minste Kvadraters Metode" eller "Oridnary Least Squares" (OLS), og går ut på å finne den linjen som er slik at summert kvadratet av den vertikale avstanden til alle prikkene, minimeres.
Vi skal ikke gå så nøye inn på matematikken her, men poenget er altså å finne en linje som passer så godt som mulig til dataene.
La oss først sørge for at vi har dataene organisert riktig. Vi må først definere y-variablen (vertikal akse) og x-variablene (horisontal akse). Det er ikke sikkert at linjen vi skal tegne bør krysse akkurat i "origo", der x og y-aksene krysser hverandre. Det er tvert i mot lite sannsynlig. Vi må derfor ha med et konstantledd, som angir hvor på y-aksen linjen skal krysse.
Vi legger til et konstantledd ved å definere en variabel 'intercept' lik 1. Linjen vi skal forsøke å finne er gitt ved
$\alpha+\beta \cdot x$
Vi skal altså forsøke å finne $\alpha$ og $\beta$ som passer best til de observerte kombinasjonene av x og y. y er forventet antall leveår:
y=g['life_exp']
pd.DataFrame(y)
... og x er BNP per capita og konstantleddet som alltid er 1:
x=pd.DataFrame(np.log(g['gdp_cap']))
x['intercept']=1
x
from statsmodels.regression.linear_model import OLS
res=OLS(y,x).fit()
print(res.summary())
Over ser vi resultatet av regresjonen. Det vi bør legge spesielt merke til er kolonnen under coef, som viser estimatene. Vi ser at $\alpha=4.9496$ og $\beta=7.2028$. Disse tallene er lagret i res.params:
res.params
Vi kan derfor bruke res.params til å plotte linjen $\alpha+\beta \cdot x$
x=np.linspace(min(np.log(g['gdp_cap'])), max(np.log(g['gdp_cap'])), 100)
regression_line=res.params['intercept']+res.params['gdp_cap']*x
ax.plot(x, regression_line,color='red')
fig
Vi ser at regresjonslinjen passer svært godt til dataene. Siden vi har brukt OLS vet vi at dette er den linjen som har minst avstand i gjennomsnitt til alle punktene. "Avstand" måles som kvadrate av den vertikale differansen. Matematisk kan vi skrive det slik:
$y=\alpha+\beta \cdot x+\epsilon$
der $\alpha+\beta \cdot x$ er linjen og $\epsilon$ er den vertikale avstanden til linjen.
Det viktigste i resulatet fra regresjonen i Eksempel 6 er ikke nødvendigvis koefisientene, selv om de er nødvendige for å tegne en regresionslinje. Enda viktigere er det å vite om helningen like gjerne kunne vært null. I så fall kan vi ikke konkludere med at det er noen sammenheng med variablene. Hovedpoenget med å kjøre en regresjon er jo vanligvis å undersøke om det er en linær sammenheng. I eksemplet over ønsker vi for eksempel å vite om det er en sammenheng mellom BNP per innbygger og antall leveår.
Når en linær sammenheng ikke kan skyldes tilfeldighet kaller vi sammenhengen for signifikant.
På grunn av ren tilfeldighet er sjeldent helningen i en regresjon null. Vi får derfor ikke noe ut av å se på selve koefisienten om den er null eller ikke. I stedet må vi se på p-verdien. p-verdien angir sannsynligheten for at resultatet skyldes tilfeldighet. Når denne sannsynligheten er mindre enn én prosent, konkluderer vi vanligvis med at resultatet er signifikant.
Koden under generer et datasett som du kan analysere med statsmodels. Kjør en regresjon med dataene, og forsøk med ulike verdier for betafor å finne ut hvor liten den må være for at resultatet ikke er signifikant:
beta=1
N=100
x=5+np.random.rand(N)
data=pd.DataFrame({'y':-3+beta*x+np.random.rand(N), 'x':x})
data
Finn en tabell på nettet med tall. Få tabellen inn i python ved å skrape nettsiden, kjør en reggresjon og kommenter resultatet